%% 函数或者脚本说明
%{
---------------------------------------------------
*文件名: ProcessHTL Version 11.19
*功   能:尝试自己处理HTL的介电参数
*注意事项:
		%

---------------------------------------------------
%}

clc; clear; clf; close all;
Input_Diretory = ".\";
FullName = @(x) Input_Diretory + x + ".csv";

%{1

Refe1 = readCST(FullName('ref01-SampleTimeDomain'), 4, 2);
Refe2 = readCST(FullName('ref02-SampleTimeDomain'), 4, 2);
Refe3 = readCST(FullName('ref03-SampleTimeDomain'), 4, 2);
Refe4 = readCST(FullName('ref04-SampleTimeDomain'), 4, 2);
% figure(7)
% plot(Refe1(:, 1), Refe1(:, 2), 'b', 'Linewidth', 1.5); hold on;
% plot(Refe2(:, 1), Refe2(:, 2), 'b', 'Linewidth', 1.5);
% plot(Refe3(:, 1), Refe3(:, 2), 'g', 'Linewidth', 1.5);
% plot(Refe4(:, 1), Refe4(:, 2), 'r', 'Linewidth', 1.5);
% legend('R1', 'R2', 'R3', 'R4', 'Location', 'northeast');

Samp1 = readCST(FullName('sam1-1-SampleTimeDomain'), 4, 2);
Samp2 = readCST(FullName('sam1-2-SampleTimeDomain'), 4, 2);
Samp3 = readCST(FullName('sam2-1-SampleTimeDomain'), 4, 2);
Samp4 = readCST(FullName('sam2-2-SampleTimeDomain'), 4, 2);
Samp5 = readCST(FullName('sam3-SampleTimeDomain'), 4, 2);
Samp6 = readCST(FullName('sam4-SampleTimeDomain'), 4, 2);
% figure(8)
% plot(Samp1(:, 1), Samp1(:, 2), 'b', 'Linewidth', 1.5); hold on;
% plot(Samp2(:, 1), Samp2(:, 2), 'b', 'Linewidth', 1.5);
% plot(Samp3(:, 1), Samp3(:, 2), 'g', 'Linewidth', 1.5);
% plot(Samp4(:, 1), Samp4(:, 2), 'r', 'Linewidth', 1.5);
% legend('R1', 'R2', 'R3', 'R4', 'Location', 'northeast');
%}

% w= {'0.823', '0.597', '0.805', '0.692', '1.689', '1.678'};

Refe = Refe2;
% Refe(:, 2) = smoothdata(Refe(:, 2), 'movmean', 15);
Samp = Samp2;
% Samp(:, 2) = smoothdata(Samp(:, 2), 'movmean', 15);
figure(10)
plot(Refe(:, 1), Refe(:, 2), 'b', 'Linewidth', 1.5); hold on;
plot(Samp(:, 1), Samp(:, 2), 'k', 'Linewidth', 1.5);
legend('Reference', 'Sample', 'Location', 'northeast');
set(get(gca, 'Title'), 'String', '时域波形');

% fft处理
FFT_refe = fft(Refe(:, 2));
FFT_samp = fft(Samp(:, 2));
Fs = 1 / (Refe(2, 1) - Refe(1, 1)); % 采样频率 THz

% T = Refe(2:end, 1) - Refe(1:end - 1, 1);
% T = Samp(2:end, 1) - Samp(1:end - 1, 1);
Length = size(FFT_refe, 1);
Frequency = Fs * (0:(Length / 2)) / Length;

figure(11)
plot(Frequency, abs(GetFFTAmplitude(FFT_refe)), 'b', 'Linewidth', 1.5); hold on;
plot(Frequency, abs(GetFFTAmplitude(FFT_samp)), 'k', 'Linewidth', 1.5);
legend('Reference', 'Sample', 'Location', 'northeast');
set(get(gca, 'Title'), 'String', '频域波形');

Smooth_r = smoothdata(GetFFTAmplitude(FFT_refe), 'movmean', 15);
Smooth_s = smoothdata(GetFFTAmplitude(FFT_samp), 'movmean', 15);
S_Final = Smooth_r ./ Smooth_s; % 相比
% plot(Frequency, abs(S_Final)); xlim([0 5]);
Phi = unwrap(angle(S_Final));
Rho = abs(S_Final);
wd = Frequency' * 1e12 * 2 * pi * 0.823e-3; % 频率X厚度
N_s = Phi * 3e8 ./ (wd) + 1;
K_s = log(4 .* N_s ./ (Rho .* (N_s + 1).^2)) * 3e8 ./ wd;

% epsilon_real = real((N_s + 1j .* K_s).^2);
% epsilon_imag = imag((N_s + 1j .* K_s).^2);
epsilon_real = N_s.^2 - K_s.^2;
epsilon_imag = 2 * N_s .* K_s;

figure(1)
plot(Frequency, epsilon_real); xlim([0.2 3]); title('real');
figure(2)
plot(Frequency, epsilon_imag); xlim([0.2 3]); title('imag');
figure(3)
plot(Frequency, angle(epsilon_real +1j * epsilon_imag)); xlim([0.2 3]); title('\delta');

function Amplitude = GetFFTAmplitude(input)
    % 取得单边谱的幅值
    Length = length(input);
    input = input / (Length);
    input = input(1:Length / 2 + 1);
    input(2:end - 1) = 2 * input(2:end - 1);
    Amplitude = input;

end
